clear all;
close all

%% Code for the paper "Existence of the Wealth–Consumption Ratio in Asset Pricing Models with Recursive Preferences"
% by Walter Pohl, Karl Schmedders and Ole Wilms

% Codes written by Ole Wilms
% Contact: www.olewilms.com

% The code numerically computes the existence condition Omega(theta) given
% in equation (18) when the states follow discrete time affine processes
% as in equation D.1 of the online appendix

% Inputs that are required are:
% - the vector of means: MU
% - the matrix F
% - the matrices h and H that determine G
% - the vector l0 and matrix l1 for the jump intensities
% - the vector of MGFs of J

% Below we provide sample code for i) the model of Bansal & Yaron (2004) as considered in Section 1.4.2
% and ii) the model of Drechsler and Yaron (2011) as considered in Appendix C.2


%% i) Example Bansal & Yaron (2004)

% Define model parameters
beta = .998;
psi = 1.5;
gamma = 10;
mu = .0015;
rho = .987;
phie = .044;
sigma_bar = 0.0078;

alpha = 1-1/psi;
theta = (1-gamma)./(1-1/psi);

%% Set up vectors and matrices for state dynamics as in Section D.1

S = 2; % Number of processes (consumption + x)

N = S+1; 

MU = [1;mu;0]; % First row is constant term, second row consumption growth, then other states

F = [0 0 0;
    0 0 1;
    0 0 rho];

h = diag([0 sigma_bar^2 (phie*sigma_bar)^2]);

H = zeros(N,N,N);
H_l = diag([0 0 0]);
H(:,:,3) = H_l;

l0 = zeros(N,1);
l1 = [0 0 0;
    0 0 0;
    0 0 0];

% Define moment generating functions for jump terms
MGF_J = @(v) [0 0 0]; 

tau = zeros(1,N);
tau(1) = log(beta);
tau(2) = alpha;

%% ii) Example Drechsler and Yaron (2011)
% 
% % Define model parameters
% 
% beta = .999;
% psi = 2;
% gamma = 5;
% 
% mu = .0016;
% phi_c = 0.0066;
% w_c = 0.5;
% 
% rho_x = .976;
% phi_x = .032*phi_c;
% l1_x = 0.8/12;
% mu_x = 3.654*phi_x;
% v_x = 1;
% 
% Es_bar = 1;
% rho_s_bar = 0.985;
% phi_s_bar = 0.1;
% 
% Es = 1;
% rho_s = 0.87;
% phi_s = 0.35;
% l1_s = 0.8/12;
% mu_s = 2.55;
% v_s = 1;
% 
% alpha = 1-1/psi;
% theta = (1-gamma)./(1-1/psi);
% 
% %% Set up vectors and matrices for state dynamics D.1
% 
% S = 4; % Number of processes
% N = S+1;
% 
% F = [0 0 0 0 0;
%     0 0 1 0 0;
%     0 0 rho_x 0 0;
%     0 0 0 rho_s_bar 0;
%     0 0 0 (1-rho_s) rho_s];
% 
% MU = (eye(N,N) -F)* [1;mu;0;Es_bar;Es];
% 
% h = diag([0 phi_c^2*(1-w_c)^2 0 phi_s_bar^2 0]);
% 
% H = zeros(N,N,N);
% H_s = diag([0 phi_c^2*w_c^2 phi_x^2 0 phi_s^2]);
% H(:,:,5) = H_s;
% 
% l0 = zeros(N,1);
% l1 = [0 0 0 0 0;
%     0 0 0 0 0;
%     0 0 0 0 l1_x;
%     0 0 0 0 0;
%     0 0 0 0 l1_s];
% 
% % Define moment generating functions for jump terms
% MGF_x = @(v_i) (1+mu_x/v_x*v_i).^(-v_x)*exp(v_i*mu_x);
% MGF_s = @(v_i) (1-mu_s/v_s*v_i).^(-v_s);
% 
% MGF_J = @(v) [0, 0, MGF_x(v(3)), 0, MGF_s(v(5))];
% 
% tau = zeros(1,N);
% tau(1) = log(beta);
% tau(2) = alpha;

%% Copmute existence condition (18)

I = 1000000; % Iteration limit

% A is the matrix for theta = 1 to compute the Jensen bound, 
% A_mink is the matrix for theta = (1-gamma)./(1-1/psi) to compute the minkowski bound 

A = zeros(I,N);
A_mink = zeros(I,N);

for i = 2:I
    A(i,:) = h_func(A(i-1,:) + tau,MU,F,h,H,l0,l1,MGF_J);
    A_mink(i,:) = h_func((A_mink(i-1,:) + tau)*theta,MU,F,h,H,l0,l1,MGF_J)./theta;
end

A = A';
A_mink = A_mink';

% Limiting behavior:

limdA = A(:,end)-A(:,end-1);
limdA_mink = A_mink(:,end)-A_mink(:,end-1);

% Existence Conditions:
ExCond_jens = (sum(limdA)) % Condition for theta = 1
ExCond_mink = (sum(limdA_mink)) % Condition for theta = (1-gamma)./(1-1/psi)


if theta < 1
    Conv= ['As ExCond_jens < 0, a model solution exists for any gamma >=',num2str(1/psi)];
    NoConv = ['As ExCond_mink > 0, there does not exist a solution'];
    Incon= ['As neither ExCond_jens < 0 nor ExCond_mink > 0, existence results are inconclusive'];

    if ExCond_jens < 0
        disp(Conv)
    elseif ExCond_mink > 0
        disp(NoConv)
    else
        disp(Incon)
    end

else
    Conv= ['As ExCond_mink < 0, a solution exists'];
    NoConv = ['As ExCond_jens > 0, there does not exist a solution'];
    Incon= ['As neither ExCond_mink < 0 nor ExCond_jens > 0, existence results are inconclusive'];

    if ExCond_mink < 0
        disp(Conv)
    elseif ExCond_jens > 0
        disp(NoConv)
    else
        disp(Incon)
    end

end
